#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "head.h"
#include "create_mesh.h"


/**
 * 将节点指标从三指标换算为单指标
 * @param i Index in x-direction
 * @param j Index in y-direction
 * @param k Index in z-direction
 * @param N Number of Mesh
 *
 * @return 节点(i,j,k)的单指标表示
 */

// 任何边界条件，都需要这种变换，因为方便并行。
int Locate_Function(int i, int j, int k, int N)
{
	return (i+2*N+1)*pow(4*N+1,2) + (j+2*N+1)*(4*N+1) + (k+2*N+1);    // 4n+1 ???
}


/**
 * @file   create_mesh.cpp
 * @author W. P. Yao <ywp@pku.edu.cn>
 * @date   2015年11月11日22:10:39
 * @brief  生成空间网格并以文本形式输出
 */

/**
 * 网格生成函数
 * @param  N 网格参数，与空间步长Delta_h满足N * Delta_h = L
 * @return   返回0表示网格生成成功，否则表示失败
 */


int Create_Mesh(int N)
{
	// 初始化fp的文件指针为NULL
	FILE *fp1 = NULL;
	FILE *fp2 = NULL;
	
	// 打开存放网格信息的文本文件
	fp1 = fopen("../Data/Point_Table.txt","w+");
	fp2 = fopen("../Data/Connection_Table.txt","w+");
	
	// 若打开失败，则网格生成失败
	if(fp1 == NULL || fp2 == NULL) return -1;
	
	// 立方体网格点
	int Domain_D = pow(4*N+3,3); 
	// 动态内存分配
	int *Table_D = (int*)malloc(Domain_D*sizeof(int));
	
	// 初始化格点，立方体D内；
	int n = 0;
	//for(int i=-2*N-1;i<=2*N+1;i++)
	//	for(int j=-2*N-1;j<=2*N+1;j++)
	//		for(int k=-2*N-1;k<=2*N+1;k++)
	for(int i=-2*N;i<=2*N;i++)
		for(int j=-2*N;j<=2*N;j++)
			for(int k=-2*N;k<=2*N;k++)
				Table_D[n++] = -1;
				printf("\n\tn = %d\n", n);
	
	// Points_Count 网格节点总数
	// int AN_Max = (3 + 16*N + 24*N*N + 32*N*N*N)/3; // 正八面体内，要改
	int AN_Max = pow(4*N+1,3);
	
	FILE *fp0 = NULL;
	fp0 = fopen("../Data/Mesh_Info.txt","w+");
	fprintf(fp0,"%d\n%d\n", N, AN_Max);
	if(fp0 == NULL) return -2;
	fclose(fp0);
	fp0 = NULL;
	
	// Create Point Table
	int AN = 0;
	
	// 以外部大正八面体为边缘；
	//for(int i=-2*N;i<=2*N;i++)
		//for(int j=-(2*N-abs(i));j<=2*N-abs(i);j++)
			//for(int k=-(2*N-abs(i)-abs(j));k<=2*N-abs(i)-abs(j);k++){
	// 以外部大立方体为边缘；
	for(int i=-2*N;i<=2*N;i++)
		for(int j=-2*N;j<=2*N;j++)
			for(int k=-2*N;k<=2*N;k++){
				
				int Domain = -1;
				int Norm1  = abs(i) + abs(j) + abs(k);

				if(Norm1 == N)
					Domain = INNER;    // 内部正八面体上的点
				else if(abs(i) == 2*N || abs(j) == 2*N || abs(k) == 2*N)
					Domain = OUTER;    // 外部立方体上的点
				else
					Domain = NORM;	   // 除此之外其他正常的点
									   // 不会有点在外部立方体之外 ???


				int Locate = Locate_Function(i,j,k,N);
											  // Index in Domain D
				Table_D[Locate] = AN;    // Define Table D and normalize.

				fprintf(fp1,"%d %d %d %d %d\n", AN++, i, j, k, Domain);
											  // Write point table
				}

	printf("\n\tAN = %d AN_Max = %d.\n", AN, AN_Max); 
	if(AN != AN_Max) {
		printf("\n\tCreate Points Error.\n");
		printf("\n\tAN = %d AN_Max = %d.\n", AN, AN_Max); 
		printf("\nEnter to continue...\n");
		getchar();
		}                      // Check


	// Create Connection Table 
	AN = 0;
	//for(int i=-2*N;i<=2*N;i++)
		//for(int j=-(2*N-abs(i));j<=2*N-abs(i);j++)
			//for(int k=-(2*N-abs(i)-abs(j));k<=2*N-abs(i)-abs(j);k++){
	for(int i=-2*N;i<=2*N;i++)
		for(int j=-2*N;j<=2*N;j++)
			for(int k=-2*N;k<=2*N;k++){
				fprintf(fp2,"%d %d %d %d %d %d %d\n", AN++,
					Table_D[Locate_Function(i-1,j,k,N)], 
					Table_D[Locate_Function(i+1,j,k,N)], 
					Table_D[Locate_Function(i,j-1,k,N)], 
					Table_D[Locate_Function(i,j+1,k,N)], 
					Table_D[Locate_Function(i,j,k-1,N)], 
					Table_D[Locate_Function(i,j,k+1,N)]);
				} // -1 is the unused point

	fclose(fp1);
	fclose(fp2);
	fp1 = NULL;
	fp2 = NULL;

	free(Table_D);
	
	return 0;

}


